Straight cracks in dynamic brittle fracture 
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We study the dynamics of cracks in brittle materials when the velocity of the crack is comparable 
to the sound velocity by means of lattice simulations. Inertial and damped dynamics are analyzed. 
It is shown that dissipation strongly influences the shape of the crack. While inertial cracks are 
highly unstable, dissipation can stabilize straight cracks. Our results can help to explain recent 
experiments on PMMA. 



I. INTRODUCTION. 

The dynamics of cracks in brittle materials such as 
glasses has recently attracted a great deal of interest. 
While an extensive body of work exists on the properties 
of quasistatic cracks, crack propagation when the crack 
grows at velocities comparable to the sound velocity is 
still poorly understood (seetl). Particular attention has 
been devoted to the study of crack_tip instabilities such as 
crack branching and oscillatiorcTEl. Typically, the crack 
tip reaches a critical velocity of the order of the Rayleigh 
speed in the material; faster cracks branch or oscillate. 
Interesting patterns were also observed under an applied 
thermal gradient!! 

In the following, we analyze crack tip instabilities in 
brittle materials. In these systems, the stress distribu- 
tion around the crack is assumed to be well described by 
the continuum theory of elasticity^. We assume that the 
instabilities observed in the experiments cited above are 
determined by these stresses. 

The stress distribution near the tip. of a moving crack 
was analytically calculated by YoffeEI. The calculation 
shows a bifurcation at a critical velocity cy sa O.Gcr 
where cr is the Rayleigh velocity. Beyond this velocity, 
the stress component tangent to the crack tip (agg(r, 9), 
assuming that the tip has radius of curvature, r) has 
a maximum at a finite angle with the crack direction. 
This result can be interpreted as a tendency for the crack 
tip to deviate from the straight direction. 

This criterion is the simplest which predicts an insta- 
bility of an inertial crack. Alternative criteria can be ob- 
tained by slightly perturbing the crack shape, and looking 
for the growth of the perturbation. InE3, a wavy pertur- 
bation is added to a straight, quasistatic crack, and the 
induced modifications of the stresses at the tip are calcu- 
lated. An instability is identified when the shear stresses 
are enough to deviate the crack from its initial direc- 
tions. In this quasistatic finite external shear is 
required to induce the instability. This analysis is ex- 
tended to dynamical cracks ina, where it is shown that 
above a certain velocity an infinitesimal shear distortion 
is amplified. The critical velocity depends on material 



parameters which describe the forces at the crack tip. 
Quasistatic cracks in PMMA under different stress dis- 
tributions follow trajectories well described-jiti terms of 
the stress intensity factor near the crack tipl!3. 

We now concentrate on cracks in PMMA, which is a 
glassy polymer. The microscopic aspects of the fracture 
process are not understood. It is possible to define a 
characteristic length in terms of the tensile strength, / t , 
the crack surface energy, Gp, and the Young's modulus 
of the material. The standard definition for thick plates 
and planar deformations ist 3 ! l ch = EGp/((l — ^ 2 )ff)- 
This length, derived from macroscopic parameters, gives 
an estimate of the scale at which continuum elasticity 
may cease to be valid. Using, the parameters in Table I, 
we find: l c h ~ 60.2 micronst 3 ! Typical structures at the 
crack tip, such as its radius of curvature, have dimen- 
sions comparable to this length. Thus, a comprehensive 
model of the growing crack should, at least, take into 
account physical phenomena beyond this scale. Elastic 
waves of wavelengths of micrometers have frequencies in 
the gigahertz range. 



c e (cal g- 1 K ) 

k (cal cm" 1 s" 1 

P (g cm" 3 ) 

l/(npc e ) 

E (GPa) 

v (Poisson ratio) 

G F (N m" 1 ) 

ft (MPa)^ 

c R ( m s 1 ) 



0.28 
4.7xl0~ 3 

1.2 
6.3xl0 3 

2.9 
0.401 

290 

130 

989 



TABLE I. Experimental values for some constants of 
PMMA relevant to the present work. The meaning of the 
symbols is: c e , specific heat, k, thermal conductivity, p, den- 
sity, E, Young's modulus, Gf, crack surface energy, ft, tensile 
strength, and cr, Rayleigh velocity. 
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In the following we will analyze cracks in PMMA by 
means of an approximate model in which details at the 
atomic and molecular scale are neglected. We take it as a 
coarse-grained approximation to a more microscopic de- 
scription in order to gain information on the role of effec- 
tive macroscopic parameters on the shape of the growing 
crack. The model has already been used in studying the 
influence of thermal gradients on crack growtlo. As we 
will discuss, we find that the viscosity (which determines 
the sound attenuation, for instance) plays a major role 
in stabilizing straight cracks and controlling their insta- 
bilities. 

In the following section, we give a brief discussion of 
the experimental situation. Then, we discuss the model, 
its general features, and related results available in the 
literature. We next present our results. The following 
section discusses possible improvements of the model, 
and the article ends with a conclusion section. There 
is an appendix where possible mechanisms which may 
lead to dissipation in PMMA are explored. 

II. EXPERIMENTAL FACTS 

The propagation of cracks, under, mode I conditions, 
in PMMA shows different regimesaQ. Most experiments 
are done in PMMA sheets under uniaxial stress. An ini- 
tial straight crack grows with a velocity dependent on 
the applied stress. Above a certain threshold, the ve- 
locity shows oscillations in time, although the crack re- 
mains straight. At higher values of the average veloci- 
ties the crack surface becomes rough due to the forma- 
tion of microscopic side branches (microbranching transi- 
tion). At even higher velocities, the crack branches into 
several paths. The transition from velocity oscillations 
(and acoustic emission) to microbranching is accompa- 
nied by a discontinuity in the average velocity. These 
transitions take place at velocities which are a fraction 
of the Rayleigh velocity, cr = 989 m/s. Only a small 
fraction of the energy dissipated during the growth pro- 
cess is radiated into sound-, wavesB. Significant heating 
effects have been reportedtj. When the growing crack is 
perturbed by means of external sound sources, many fea- 
tures of the previous picture are modified!!]. The velocity 
gap at the microbranching transition is washed out. 

The previous, picture also seems to hold for cracks in 
ordinary glasal3. Cracks moving at constant speed seem, 
however, much more difficult to stabilize in ordinary glass 
than in PMMALLia. 



III. DISSIPATION IN PMMA 

Elastic waves are attenuated in real materials, and en- 
ergy is transferred to degrees of freedom other than those 
which describe sound waves. This attenuation can be 



modeled by adding a viscosity term u to the elastic equa- 
tions of motion of the form r/W 2 dtU where u is the dis- 
placement, and 77 is a viscosity coefficient. In this long- 
wavelength limit transverse sound waves acquire an at- 
tenuation a = rjk 2 /2pcT where p is the density and ct is 
the transverse sound velocity. Thus there is a wave- vector 
at which the attenuation of a wave becomes comparable 
to its wavelength, a A = irrjk/ pcx ~ 1. Beyond this scale, 
sound waves are overdamped, and the analysis reported 
irJ3 certainly needs to be modified. 

The influence of a different form of viscosity on the 
velocity of straight cracks was considered inlia. It was 
found that the presence of damping at the edges of a 
type III crack leads to a steady state velocity which, at 
high viscosities, is inversely proportional to the damping 
coefficient. 

The term 77 V 2 dtu is thought to be appropriate at very 
low frequencies. However, in glassy systems the attenua- 
tion is a complicated function of frequency due-.to the 
different relaxation processes which contributes. For 
example, iat-PMMA at high frequencies (several GHz), 
a A ~ 0.lEic3, and, thus, a cc uj. At lower frequencies 
(2 MHz), the dependence of a on frequency can be fit- 
ted by a power law, a ~ w c , with c ~ 1.5 — 1.7c3. It 
has been argued that some relaxation processes freeze 
below 165KO. It is likely that, at frequencies < 100 
GHz a dependence other than uj 2 may arise. At suffi- 
ciently low temperatures, the situation simplifies some- 
what, as the main modffi which contribute to dissipation 
are better understood^. A microscopic analysis of dis- 
sipation processes at low temperatures, using as input 
detailed experimental data on the low energy modes in 
glassy polymerst£l, is given in the Appendix. 

IV. GENERAL FEATURES OF CRACKS IN 
BRITTLE SOLIDS 

A. The Equations of Motion 

The equations of motion including viscous terms for an 
isotropic medium can be written as: 

pd tt u = MV 2 u + (A + /Lt)V (Vu) + 

77V 2 a t u + (V> + t?)V (V0 t u) , (1) 

where u is the displacement field, A and /1 the Lame 
coefficients and ip and 77 the corresponding coefficients in 
the viscous case. 



B. Attenuation of Elastic Waves 

The relation between the attenuation coefficient and 
the frequency can be derived approximating the solutions 
of the equations of motion by a longitudinal plane wave, 

u(r,t) = uoe i[(k+iot)x ~' ,rtl (2) 
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where Uo = uqi is the amplitude of the wave and i the unit 
vector in the x-direction. The result for the attenuation 
coefficient is, 



y/(X + 2 M ) 2 + (i/j + 2f]) 2 U 2 
(A + 2[i) 2 + (-0 + 2t]) 2 uj 2 



(3) 




FIG. 1. Attenuation coefficient divided by the frequency 
as function of the frequency obtained for longitudinal plane 
waves (see text). The frequency is expressed in units of the 
frequency at which the maximum occurs (a> ma x). 

The dispersion relation (k(u>)) is just Eq. ^| changing 
in the right hand side the minus sign by a plus sign. 
In Fig. [y we plot a/cu versus uj/uj maK . This ratio has a 
maximum at, 



= V3 



A + 2/i 



(4) 



ip + 2?] 

and the behavior of the attenuation coefficient in the low 
and high frequency limits is, 



u) — > 0, a oc lu 



oo, a oc oj 



1/2 



(5a) 



(5b) 



These results show that, roughly, the attenuation coeffi- 
cient saturates at wavelengths at which the attenuation 
and the frequency of the wave become comparable. On 
the other hand, the behavior at low frequencies is con- 
sistent with the experimental data and with the results 
obtained by means of the microscopic analysis described 
in the Appendix. At high frequencies, however, both the 
experiments and the microscopic analysis give a oc u>. 
This discrepancy is not surprising, considering that the 
continuum theory should fail at small length scales (see 
above). 



C. Generalization of the Griffith criterion 

The Griffith criterion is a fundamental element in the 
theory of fracturecll. According to Griffith, a crack starts 
to advance if, in increasing its length by 6L, the elas- 
tic energy released is greater than the amount of energy 
needed to create the new fracture surface. Mott wa 
first to include kinetic effects in Griffith's analysis 
He proposed to add a kinetic energy to the Griffith's to- 
tal energy. However some of the conclusions inferred from 
Mott's analysis are not valid, for instance the predicted 
value for the maximum crack speed was lower than ex- 
pected (seed for a full discussion). 

Here use a different approach and attempt to directly 
generalize the Griffith criterion to the viscous case by 
keeping track of the energy flows. We balance the energy 
release, the difference of the elastic energy and the sur- 
face energy, with losses due to viscous dissipation. We 
consider a long system of width W and thickness d, with 
a crack of length L. To estimate the elastic energy release 
E r , we note that for L << W a roughly round region of 
diameter L is fully relaxed, so that E r oc e 2 L 2 d where 
e is the strain that causes the material to break. For 
L >> W we must put E r ~ e 2 WLd. The second term, 
the cost of creating new fracture surface Ef, is Ef oc Ld. 
Finally, the dependence of the rate of viscous dissipation 
5 Ed on the crack speed v, may be estimated for slow 
cracks from a symmetry argument: Since 5 Ed — > as 
v — > 0, and must be non-negative for any v, we conclude 
that 5 Ed oc v 2 St. The coefficient of this term goes to 
zero as r\ — > 0, so that we put Ed oc rjv 2 St. Now we set 
SEd = 5E r — SEf, and use 8L — vSt. For short cracks we 
see that 



r\v oc e L — q 



(6) 



where q is a constant. Short cracks accelerate. For long 
cracks, on the other hand, there is a terminal velocity: 



rjv oc e W 



(7) 



If the terminal velocity is less than cy, the Yoffe thresh- 
old, we may expect that the crack will never branch. 

Note that the analysis of this section can be expected 
to be valid only in the limit of low crack speed. In partic- 
ular, the viscous dissipation term on the right hand side 
of Eq. involves the motion of the lattice in response to 
a passing crack. Thus the rj in this equation is not quite 
the same as the one above, and for substantial v proba- 
bly has a complicated dependence on the microscopic rj 
and on v itself. We will test these ideas with simulations, 
below, and find that for speeds << cr Eq. [?] is rather 
well obeyed, but that there are deviations at large speed. 



D. The branching instability 

The assumption in the previous section is that the crit- 
ical speed for branching is independent of rj, which is 
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what we find numerically for small rj (see below). This 
is a bit unexpected since, in the presence of dissipation 
the analytical solution of Yoffe, for exampleja is no longer 
correct. We can see where this assumption would break 
down by examining the form of that solution. 

The stress field described can be derived from an ap- 
propriate distribution of forces applied at the crack edges 
which have the general form f(r — vt). The stresses 
at an arbitrary point of the plane can be obtained by 
means of the Green function, G{j(r — r',t — t'), with 
Fourier transform Gy(k, lu). In the absence of dissipa- 
tion, the frequency u) appears only in combinations of 
the type J 2 — c 2 LT k 2 , where cl,t denote the velocity 
of longitudinal and transverse sound waves. Dissipation 
changes these expressions into lo 2 — c 2 L T k 2 — ir/ojk 2 / p. 
The Fourier transform of the applied forces can be writ- 
ten as f(k, lu — v • k). Hence, the denominators in the 
Green's functions become (v • k) 2 — c 2 L T k 2 — ir/k 2 v ■ k/p. 
At low values of fc, the influence of the viscosity is negli- 
gible. The long distance behavior is well described by the 
Yoffe solution. For large k on the other hand, the viscous 
term dominates. This term is more anisotropic than the 
inertial term, as it contains one power of v • k, instead 
of two. Hence, we expect the tendency towards insta- 
bility to be reduced. We can see when this is relevant 
by putting k ~ 1/a, and noting that large k means that 
k pv/rj, or equivalently, aa >> 1. For very large 77 the 
branching threshold should eventually shift. As we will 
see, we have confirmed this shift using the simulations. 



V. NUMERICAL SIMULATIONS 



A. Discretization of the Equations of Motion 



the coupling remains zero at all latter times. The model 
used here is deterministic, and the system is always out 
of equilibrium. 



B. The elastic constants of the model 

In order to find the relationship between the parame- 
ters of our model and the experimental constants we need 
to write down the equation of motion of our model (Eq. 
(B)) in the continuum limit, 



mduVL = Ka 2 



^V 2 u+^V(Vu) 



770 a 



§v 2 9 iU + ^v(va t u) 

o 4 



(9) 



where a is the lattice spacing. Solving the above equa- 
tions in a finite difference scheme gives Eq. (||). Then, 
comparing Eqs. ([!]) and (j^) one can obtain the relations 
we are seeking for. First note that our model gives the 
following relations between the continuum parameters, 
A = /i = 3-E/8 and ip — rj and a Poisson coefficient of 
1/3. On the other hand, 



3K 

8m 



3E 



(10) 



where a is the lattice constant of our triangular spring 
network. The longitudinal and transverse sound veloci- 
ties are, 



cl 



Vs. 



ct 



' A + 2/i 



P 



(11) 



In our simulations we work in two dimensions, and 
discretize the continuum equations of elasticity by us- 
ing a spring model on a triangular two dimensional ja 
following our previous work on quasistatic fractureo 
The equation of motion for the displacement, u r of the 
node at r combines inertial and viscous terms. In our 
discrete model, we get the k 2 dependence of the attenua- 
tion by using the fact that the friction forces can depend 
only on the relative velocities of neighboring nodesa. The 
equations of motion are: 



From these results the Rayleigh speed can be easily 
derived!! 



c R = 0.9325c T 



(12) 



Finally, the equations which relate the constants of our 
model with the macroscopic parameters of the material 
are: 



m = 3pa 2 d/8 



(13a) 



m - 



d 2 u r 

dt 2 



= ^Kh[h- (u r - u r /)] 



du r > 



(8) 



where the sums in the second term are over the nearest 
neighbor nodes, r', to node r and n is the unit vector 
from node r', to node r . The fracture process is de- 
scribed by deleting the forces between two nearest neigh- 
bor nodes when the relative strain, |n • [u r — u r '] |, ex- 
ceeds a threshold, u t h- This process is irreversible, and 



K = 8c^m/3a 2 



Tj = 7]d 



(13b) 



(13c) 



where d is the thickness of the sample. In the next sub- 
section, we discuss some difficulties in relating uth to 
macroscopic variables, particularly due to the discretiza- 
tion scale a. 
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C. Fracture threshold and macroscopic variables 

In the following we take Utn — 0.1a. The process of 
failure is described in macroscopic terms by the maxi- 
mum load above which the material fails. In our model, 
the maximum force exerted by the springs, Ku t h, should 
be equal to this load times the lattice spacing, a, times 
the thickness of the slab, d. As K is independent of a, 
we find that Uth °c a. 

The description of fracture by discrete element meth- 
ods always leads to failure at the smallest scale. In our 
case, a sample under load will eventually fail through the 
snapping of a single row of springs, and crack widths are 
always of order aE3. The energy required to create a crack 
is, because of this effect, dependent on the discretization, 
and goes to zero as the lattice spacing is decreased. In 
the present case, the energy needed to create a crack of 
(macroscopic) length I goes as ^Ku^ h oc a. Thus, we 
cannot fit the tensile strength and the fracture energy 
at the same time in a discretization independent way. 
This drawback can be inferred from the existence of a 
characteristic length which combines these quantities, as 
analyzed in the introduction. 

Macroscopic crack energies can be obtained from mod- 
els which incorporate non local. effectSLJ. However, de- 
tailed microscopic simulationsLj show cracks of atomic 
width, with little or no damage outside a zone of mi- 
croscopic dimensions. They are also difficult to reconcile 
with the existence of crack energies with macroscopic val- 
ues. The origin of macroscopic failure zones in quasibrit- 
tle materials such as PMMA is not well understood. 



D. The dissipation term 



the level of discretization. Only for a given lattice con- 
stant, a, the crack energy and the maximum load can be 
consistent with macroscopically determined values. On 
the other hand, instability criteria based on energy con- 
siderations, such as the Griffith criterion and extensions 
thereof, depend only on energy differences. Hence, we 
do not think that the problem discussed here is a seri- 
ous obstacle to the analysis of the instabilities of moving 
cracks. As there are substantial differences between dis- 
crete and continuum modelsEj, it would be interesting 
to analyze further the relevance of the intrinsic length 
scale determined by the macroscopic parameters which 
describe fracture. 



F. Numerical Procedures 

Simulations have been performed in rectangular strips 
with the y orientation along one of the axis of the tri- 
angular lattice. The lattices shown in the figures are 50 
nodes wide and 275 nodes long. The boundary condi- 
tion is fixed displacement of the edges so that the initial 
strain is below uth- Then, bonds near the lower horizon- 
tal edge are broken at a fixed rate, so that the velocity 
of the crack is well below or- To integrate the equa- 
tions of motion we use Heun's method with a time step 
small enough for appreciate small variations in the ve- 
locity of the crack. Once the crack is sufficiently long, 
strains near its tip begin to exceed Uth, and the crack 
continues growing by itself. Very short cracks do not 
grow on their own, because the strains at the tip do not 
exceed u t h- The minimum size for self-sustained growth 
decreases with increasing dissipation. 



Qualitatively, each node represents a region in the ma- 
terial comparable to the scales relevant to the experimen- 
tal situation. In our case, we have l c h ~ 60 microns. We 
use a phenomenological damping term, r\ ~ 1 in units 
where K — m = 1 , which implies that sound waves at this 
length scale are in the overdamped regime. As mentioned 
earlier, the sound attenuation in glassy polymers has a 
complicated dependency on frequency. Our choice of 
overestimates the experimental value in the GHz range 
although probably underestimates it at lower frequencies 
(note that our model assumes that oA ~ A -1 at all wave- 
lengths). The scheme used here is intermediate between a 
full scale atomic simulations, and more phenomenolog- 
ical modelsEj, where dissipation takes place within the 
units used in the discretization. 



VI. RESULTS 

A. Stable and unstable cracks 

In the absence of damping, straight cracks become un- 
stable on short time scales. Typical results for cracks 
growing in a narrow slab under an applied strain at the 
edges are shown in Fig. || (a)-(c). The crack tips acceler- 
ate, exactly as predicted in Eq. (||) until they approach 
cy, and then they branch. The velocity of the upper- 
most part of the crack pattern is depicted in Fig. ^ (d)- 
(f). We note that the crack velocity strongly oscillates as 
a function of the crack length. 



E. Drawbacks of the model 



As mentioned earlier, the main difficulty with the 
model is the fact that the crack energy does not scale with 



5 




200 100 200 
y (lattice param.) 

FIG. 2. Behavior of inertial cracks when the external 
strain, u app i, is varied, (a) u app i=0.02'i. (b) u app ; =0.026. (c) 
Uappi— 0.028. And also their velocity (in units of or) as func- 
tion of position of the advancing crack tip. (d) u app i =0.024. 
(e) u app i=0.026. (f) it aj) j,i =0.028. The threshold for breaking 
is u t h = 0.1. 

We find that straight cracks can be stable (cf Fig. |^a, 
below) , and move at constant velocity, in the presence of 
dissipation. As the driving force is increased, we observe 
a branching instability. This behavior is what is pre- 
dicted by Eq. (fjj). If the terminal velocity is below cy 
(which we assume to be independent of rj, see below) the 
crack will be slowed down and prevented from branch- 
ing. This behavior is shown in Fig. || where the terminal 
stable velocity in units of cr is plotted against the dis- 
placement at the borders of the system. The evolution to 
increasing velocity proceeds as the external displacement 
is increased. The curves end at the branching instability. 
From them, we can deduce that this threshold is inde- 
pendent of the parameters of the simulation (size and 
viscosity) and ss 0.7cr, in agreement with Yoffeu calcu- 
lations. 

The reason for the inscnsitivity of the branching 
threshold to r\ was given above. We find in our sim- 
ulations that to shift the threshold significantly below 
0.7cr, rj must be greater than 7 for a 50 x 300 mesh, 
which is, we think, an unphysically large value. Never- 
theless, the branching threshold reported ino ~ 0.45cr, 
is around 35% smaller than in the numerical simulations. 

The fact that we do not see an abrupt change in the 
velocity before the branching threshold is related with the 
method we use to perform the simulations. Boudet and 
CilibertoES demonstrated experimentally that this jump 
was not present when sound was added into the system. 
This is what we do in the simulations: the slow cutting 
of bonds that initiate the fast failure is a source of sound 



into the system. 
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FIG. 3. Terminal velocity for a stable crack in units of the 
Rayleigh speed as function of the applied strain for different 
values of the viscosity and width of the system. 

We can directly verify the validity of Eq. (0) by consid- 
ering a number of different sets of the parameters e, rj, W, 
such as those shown in Fig. ^ and viewing our data in the 
form of a data collapse. This is done in Fig. |I] where we 
show that rjv is very accurately a linear function of e 2 W 
for low speeds. For larger speeds, of the order of cy there 
are deviations from Eq. [j], as we expect. 
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FIG. 4. Data collapse of r/ v plotted against e 2 W. 
o,rj = 0.5; +, rjo = 1-0; □, r\ = 1.5; X, r\ = 2.0, for 
12 < W < 50 and 0.029 < e < 0.058.The straight line is 
fitted to all the points for which v < 0.3. 



B. Thermal noise and the branching threshold 

Numerically is possible to obtain lower critical branch- 
ing velocities by adding a random noise in Eq. ^|, asnwill be 
discussed in detail in a forthcoming publicationE3. This 



G 



random noise is in the form that it has zero mean at its 
correlation is 



(14) 



when i = j, and zero otherwise. T is the temperature, 
and £ is a parameter that cpntrols numerically the am- 
plitude of the noise. In Ref. E3 £ is related to r\ using the 
fluctuation-dissipation theorem. Here we take a more 
phenomenological approach in order to simply illustrate 
the effect. 

Let us take one of the stable cracks whose velocity is 
plotted in Fig g. The one with AW/W=0.Q3, r)=l, and 
L=50 has a velocity of w 0.65cr. Fig. || show the shape 
of the cracks for different values of £. The temperature 
is set to a fixed arbitrary value of 100. When the noise 
is very low (10 -9 ) nothing happens to the crack. When 
it becomes higher, always at the same terminal veloc- 
ity, small side branches like in the experiments can be 
obtained. The spacing between branches decreases as 
the noise intensity is increased, and, in some cases (see 
Fig. ||(c)), the crack can disestabilize at later stages in the 
growth process, until finally the crack becomes unstable 
(Fig. |(d)). 

Since in our simulations we have all the information 
of the displacements at all points of the network, we can 
compute the stress tensor at any time of the evolution of 
the crack. In particular it is interesting to see the values 
of (Tyy — (fax around the crack tip. According to Coterell 
and RiceO, this quantity describes the stability of the 
cracks. It should be less than zero at the crack tip for 
an stable crack and greater than zero when it becomes 
unstable. This parameter is shown in the lower part of 
Fig. H when the corresponding cracks in the upper part 
of the figure have advanced half the length of the system. 
The three cracks with lower noise intensity are, according 
to this criterion, still stable, but the one with higher noise 
has a tendency to continue deviating to the left. 



(a) 



(b) 



(c) 



(d) 



III 
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FIG. 5. Shapes of the cracks when noise is added (a) 
e = 1CT 9 , (b) e = 4xl(T 8 , (c) e = KT 7 , and (d) £ = 2xl(T 7 . 
The lower figures (e-h) show the sign of the Coterell and Rice 
parameter (a yy — a xx ): dark grey positive and light grey neg- 
ative. Stresses are taken when the crack has advanced half 
the length of the system. The noise is the one of the corre- 
sponding above figure (a-d). 

This method for representing the stress tensor with 
a wide variety of parameters in an elastic medium can 
provide a valuable tool for inspecting the conditions for 
the stability, and the analytical approximations that can 
be madea. 



VII. HEATING AND ENERGY DISSIPATION 



In the previous section, examples of crack propagation 
in the presence of thermal noise due to an external envi- 
ronment, have been discussed. However, the local tem- 
perature around the crack tip must also take into account 
the energy dissipated by the viscosity term that we have 
in the equations. Near the crack tip typical deviations 
of the nodes from equilibrium are of order a. Typical 
velocities are of order Ku t h/Tj . The energy dissipated 
per node and per unit time is ~ K 2 u 2 h /i] . In terms 
of macroscopic quantities, the energy generated per unit 
time and unit volume is ~ p 2 c 4 a 2 a 2 /(iJ,E 2 ), where c is 
some average of the longitudinal and transverse sound ve- 
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locities, E is the Young's modulus and a c is the macro- 
scopic elastic limit. This dissipation generates thermal 
gradients. They will be determined by the condition: 



dT - m 1 d£ 

— = k.V 2 T+ — 

at pc e at 







(15) 



where k is the thermal diffusion coefficient, S is the en- 
ergy being dissipated and c e is the specific heat. Assum- 
ing that most of the dissipation takes place at distances 
from the crack tip comparable to its radius, we find that 
the temperature increase at the tip can be written as: 
Altjj, ~ c 4 a 4 a^/(Kfj,c e E 2 ) This expression is highly sen- 
sitive to the choice of a, the tip radius. Hence, it is 
difficult to make accurate estimates of the expected heat- 
ing. Experimentally, significant increasesJn temperature 
near the crack tip have been reportedEj. Energy dis- 
sipation has also been observed incl, where it is argued 
that most of the energy is spent in increasing the crack 
surface. However, even for slow, straight cracks, a signif- 
icant rise in energy dissipation as function of velocity is 
reported. In our simulations the elastic energy lost when 
one spring is cut goes into surface energy, whereas the 
viscous dissipation goes into heat. Heating of the crack 
tip increases thermal noise there. This could be quite sig- 
nificant since, near, but below the branching speed the 
stress distribution becomes nearly isotropic, so that rela- 
tively small thermal effects could lead to branching. The 
considerations in this section will be worked out in detail 
in Ref£3. 



VIII. CONCLUSIONS 

The analysis reported here indicates that viscous ef- 
fects change significantly the propagation, and instabili- 
ties, of cracks in brittle materials. The general features 
that we have found should be reproduced, for example, in 
PMMA, even though the viscosity is a more complicated 
function of frequency than the one considered hei:c._5pmc 
of the characteristics of the experimental resultsodlj'Q are 
already qualitatively described by the present approach. 
In particular, our approach would explain why experi- 
ments in glass are harder to perform than in PMMA: 
its associated viscosities are lower than in PMMA and 
thus it is closer to the instability. On the other hand, the 
branching threshold seems to be lower in the experiments 
than in our numerical simulations for the chosen param- 
eters. This fact has also been discussed and shown to be 
due to other effects not contained within the model, but 
that can be incorporated as external noise. Of course, 
the richness and complexity of fracture in these mate- 
rials will require further investigations. We hope that 
the approach herewith proposed will help to improve our 
understanding of these interesting phenomena. 
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X. APPENDIX 

The energy dissipated by sound waves goes into other 
excitations of the system. Taking the sound velocity of 
PMMA to be v s ~ 10 5 cm/s, a mode of wavelength of 
one micron has a frequency of 10 9 Hz. In energy units, it 
corresponds to w 7 x 10 _7 eV. At sufficiently low temper- 
atures, the main source of inelastic scattering at these 
frequencies, as seen by neutron scattering, is quantum 
tunneling between equivalent coafigurations of the CH3 
groups attached to the polymernj. The possibility that 
these excitations play a role in sound attenuation was 
suggested ir£3. 

An acoustical phonon modulates the distance between 
polymers. At the low frequencies involved, the backbone 
of the polymer cannot be excited, and can be considered 
rigid. The motion of a nearby polymer changes the po- 
tential acting on a given CH3 group, breaking the initial 
threefold symmetry. The asymmetry in the potential in- 
duces transitions between the quantum levels of the CH 3 
group, and leads to dissipation. 

The interaction between the neutral CH3 unit and 
other parts of the polymer arise from mutual induced po- 
larization. Assuming that neither part has a finite elec- 
tric dipole, the interaction energy is given by the van der 
Waals expression: 



E : 



\ " 



\d a 



3«rO(4rVr 2 | 



„2|2 



(16) 



where r is the distance, e is the dielectric constant due 
to the rest of the material, = (0\r\m) represents a 
matrix element of r between states of unit a ( and a 
corresponding expression for d\), and A m is the energy 
difference between the ground state of unit a and a given 
excited state. 

The order of magnitude of E in ( [l6|) is E ~ e t r t A , 
where d goes as the dimension of the unit, A is a typical 
electronic excitation energy, and r is the separation. 

The splitting between configurations of the CH3 unit 
goes as the difference in the interaction energy (|l6| ) at 
nearby H sites. Hence, it goes as ^Srdri-n- 

A phonon which induces displacements ilk in a given 
molecule leads to a change in the intermolecule distance 
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r of order (kr)u k . Using the golden rule, the energy per 
unit time absorbed by a given CH3 group goes as: 



dE 
~dt 



e 4 d 4 



e 2 r 6 A 



l H—H , 2 2 



k U k U p hPtunn(iOph) (17) 



where ptunn(u) is the density of tunneling centers of en- 
ergy splitting io. The energy dissipated per unit volume 
can be estimated from (|17j) by multiplying that expres- 
sion by the number of CH3 groups per molecule, N, and 
dividing by the volume of the molecule, £1. 

On the other hand, the energy dissipated per unit vol- 
ume, can be written as& 



dE , ■-, , 

— OC 7]UJ ph k u k 



(18) 



where rj stands for an average of the macroscopic viscos- 
ity. In terms of 77, the sound attenuation goes as 
where p is the density, and c is the sound velocity. Mak- 
ing use of the microscopic parameters, this gives: 



a 



fc- 



N 



d 2 

H — H 



TipVlc 2 



e 2 r 6 A 



PtunniUph) (19) 



This expression is very sensitive to the value of r, the 
distance between nearby chains. A rough estimate can 
be made by assuming ~ A ~ leV, cIh-h ~ d ~ lA, 
r ~ 5A, p(w) ~ (1/ieV) -1 and e = 1. Using these param- 
eters, we obtain a A ~ 1 — 10, where A is the wavelength of 
the phonon. The dependence on frequency of aA is that 
of Ptunn{^), which, in the relevant range of frequencies, 
~ ltf£]V, is roughly constantEj. This gives aaw, in line 
withE3. Below lower frequencies, ~ 1kHz, ptunn(u) oc u>, 



and 1 



similar to the behavior 



leading to aA oc A - 1 , 
reported inEj. 

Quantum tunneling is suppressed by thermal fluctua- 
tions, which break the degeneracy of the three potential 
minima seen by the CH3 groups. We assume that these 
fluctuations arise from changes in the van der Waals in- 
teractions with neighboring polymers. The presence of 
a lattice vibration, of momentum k and amplitude Uk 
induces a splitting between equivalent minima of: 



e 4 d 4 d H - 



H-H 



e 2 r 6 A 



ku k 



(20) 



The derivation of this expression follows the analysis pre- 
sented earlier. At finite temperatures, u k shows random 
fluctuations. The mean square deviation of AE is given 
by: 



((AEf 



e 4 d 4 



-(k 2 u\ 



(21) 



and: 



(k 2 ut) 



n / d 3 k k kBl 



k R T 



P nn 2 c 2 k 2 n P c 2 



T 

UJD 



(22) 



where loo is the Debye temperature. In terms of the 
attenuation rate calculated above, we can write: 



((AE)' 



aA ksT 

N Ptunn(E) 



T 



(23) 



Using the same parameters as above, and setting up = 
300K, we find that the fluctuations AE are of order IpeV 
when T ~ 10K. 

At higher temperatures, the dynamics of CH3 groups 
lose coherence. They still influence sound attenuation, 
because they mediate interactions between phonons. 

Note, finally, that the mechanism discussed here does 
not lead to plastic effects. The material remains brittle at 
the scales of interest, as it seems to be the case in PMMA. 
Other mechanisms, such as the irreversible relaxation of 
structural defects, may lead to viscoplastic effectsLj. 
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